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Abstract. Some dynamical properties of a bouncing ball model under the presence 
of an external force modeled by two nonlinear terms are studied. The description of 
the model is made by use of a two dimensional nonlinear measure preserving map on 
the variables velocity of the particle and time. We show that raising the straight of 
a control parameter which controls one of the nonlinearities, the positive Lyapunov 
exponent decreases in the average and suffers abrupt changes. We also show that for 
a specific range of control parameters, the model exhibits the phenomenon of Fermi 
acceleration. The explanation of both behaviours is given in terms of the shape of the 
external force and due to a discontinuity of the moving wall's velocity. 



1. Introduction 

The bouncing ball model consists of a classical particle of mass m which is confined to 
bounce between two infinitely heavy and rigid walls [1] . One of the walls is assumed to 
be fixed while the other one moves in time according to a periodic function. This model, 
also known as the Fermi-Ulam model (FUM), is a simple dynamical system that can be 
modeled using the formalism of discrete mappings. Moreover, many tools developed to 
characterise such a model show to have great applicability in more complex mappings. 
Considering the FUM, many results are known in the literature. Particularly for the 
particle suffering elastic collisions with either walls, it is known that the phase space of 
the model is of mixed kind [2] in the sense that depending on the combination of both 
the control parameters and initial conditions, invariant spanning curves limiting the 
size of chaotic seas and Kolmogorov-Arnold-Moser (KAM) islands can all be observed. 
The presence of the invariant spanning curves yields in a limit for the energy gain of a 
bouncing particle, thus the Fermi acceleration (unlimited energy gain of the particle) is 
not observed. A similar version of the model, called as bouncer [3], consists of a classical 
particle, in the presence of a constant gravitational field, suffering elastic collisions with 
a periodically moving platform. The returning mechanism of the bouncer model, a 
mechanism that injects the particle for a next collision with the moving wall, is rather 
distinct of the FUM. In the bouncer, it is due only to the gravitational field while 
in the FUM it is given by a collision with a fixed wall. These differences yield in a 
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profound consequence for the dynamics of a bouncing particle. Depending on the control 
parameter, the unlimited energy gain is observed in the bouncer model, a phenomenon 
that is not present in FUM with periodic and smooth oscillations. The differences were 
clarified by Lichtenberg, Lieberman and Cohen Hybrid versions of the FUM and 
bouncer were recently studied [3, E] as well as a stochastic version of the FUM [7j. 

There are also many important results concerning the inclusion of damping forces on 
both the models (see for example Ref. [8j for a short review). One of them is the presence 
of a drag force [9] , so that the particle is moving inside a gas with the dissipation acting 
on the particle along its trajectory. The dynamics of the problem is, generally, given by 
a nonlinear mapping that is obtained via the solution of Newton's law. A different kind 
of dissipation can be introduced via inelastic hits of the particle with the walls. Thus, 
there is a restitution coefficient that makes the particle experiences a fractional loss of 
energy upon collisions. Despite both kinds of damping often occur in nature, they have 
profound and different consequences in the dynamics of the models. As an example, in 
Refs. p!Ol [TT] and considering inelastic collisions, Tsang and Lieberman considered the 
simplified FUM (both the walls are fixed but the particle changes energy and momentum 
upon collisions with one of the wall as if the wall were moving) with inelastic impacts. 
They have evidenced contraction on the phase space and in particular, observed the 
presence of a strange attractor. Recently, a rather similar version of the dissipative model 
[T2] . confirmed the property of area contraction and in addition, a boundary crisis was 
characterised. Additionally, a family of boundary crisis was observed when collisions 
with the two walls are inelastic [I3]. The bouncer model was also considered under 
inelastic collisions. For example, in [H] Holmes discusses the appearances of horseshoes 
in the inelastic bouncer and gave an illustration of a homoclinic orbit in such a model. 
After that, Everson [15] presents and discusses with many numerical simulations the 
appearance of period doubling cascade in the damping bouncer model. Period doubling 
cascade was also observed in [16] for the completely inelastic collisions. The presence of 
frictional force however was considered by Luna-Acosta [17] and Naylor, Sanchez and 
Swift [IE] in the bouncer model. They too observed period doubling cascades and in 
special Luna-Acosta [17] has achieved analytically dimensional reduction for the limit 
of high dissipation. 

In this letter, we study a non dissipative version of a bouncing ball model seeking 
to understand and describe some of its dynamical properties considering however that 
the motion of the moving wall is given via a crank-connecting rod scheme. For such 
a scheme, it is known that there are two nonlinearities present in the model and 
each of them play important rules in the dynamics. Depending on certain ranges of 
control parameters, there can be profound consequences on the dynamics of the system. 
Particularly, when one of the two control parameters is raised, the positive Lyapunov 
exponent experiences a drastic reduction. It is also important to say that the particle is 
in the total absence of any external field. Other important result for this model is that 
it yields, for specific control parameter values, the phenomenon of Fermi acceleration 
(unlimited energy growth). The phenomenon is characterised, for the first time in the 
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present model, in terms of a discontinuity of the derivative of the wall's position with 
respect to the time, thus leading the particle to acquire unlimited energy gain. 

This paper is organised as follows. In section [2] we present the model and the 
expressions of the mapping that fully describes the dynamics of the system. Section [2] 
is also devoted to a discussion of the numerical results and the behaviour of the positive 
Lyapunov exponent. In section [3] we propose a simplified version of the model and study 
the behaviour of the average velocity as function of a control parameter. We show that 
Fermi acceleration emerges naturally from the deterministic dynamics of the model for 
specific control parameter values. Final remarks and conclusions are drawn in section 

SI 

2. The model and the mapping 

The model is described using a two dimensional mapping for the variables where 
Vn and tn are the corresponding velocity of the particle and time immediately after the 
n^^ collision with the moving wall. We assume that one wall is fix ed a.t x = I and 
that the motion of the moving wall is given by s{t) = Rcos{wt) + \J L'^ — B? sin^ {wt) 
(we stress the term wt in the equation of s{t) is represented by the variable in Fig. 
[1]), where R denotes the radii of the crank, L is the length of the connecting rod and 
w is the corresponding frequency of oscillation. Figure [1] illustrates the model under 
consideration. We stress the equation for s{t) is easily obtained from the condition 
of -Rsin(0) = Lsin('?/'), as can be seen in Fig. [TJ Before we write the equations of the 




< > 



Figure 1. Illustration of the model under consideration. 

mapping, let us first discuss the initial conditions. We assume that, at a time t = tn, the 
particle is at the position Xp{tn) = s(t„) with velocity f = f„ > 0. Thus such an initial 
condition can be considered as if the dynamics were already running in the system along 
the time. We emphasise that two different kinds of collisions can be observed namely: 
(i) multiple hits with the moving wall and (ii) a single hit with the moving wall. In case 
(i), the particle suffers a collision with the moving wall but then, before it leaves the 
collision zone, which is defined as a; G [—R,R], the particle experiences a second and 
then successive impact with the moving wall. Such kind of collisions becomes rare in the 
limit of high energy but they are quite often to be observed in the regime of low energy. 
It is also easy to see that there are too many control parameters in the model, 4 in 
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total, namely R, L, I and w and that the dynamics of the system does not depend on all 
of them. It is convenient to define dimensionless and more appropriated variables. We 
define e = R/l, r = R/L, Vn = Vn/ (wl) and measure the time in terms of the number of 
oscillations of the moving wall 0„ = wt„. For the dimensionless variables, we consider 
that the range for e is e G [0, 1] and for r is r G [0, 1]. The limit of r ^ corresponds 
to L — oo and r — > 1 is obtained for L —>■ R'^. With this new set of variables, the 
mapping that describes the dynamics of the system is written as 

bn+i = [4>n + AT„] mod(27r) 

_j_ rcos(0,i+i) ) (1) 

where the expressions of V* and AT„ depend on the kind of collision. For case (i), which 
corresponds to the multiple hits with the moving wall, the corresponding expressions 
are V* = —Vn and AT„ = 0c with (pc obtained by the solution of G{(f)c) = with G{(j)^ 
given by 



T : 



K+i = K:-2esin 



G{(f)c) = e cos(</)„ + (f)c) -e cos(0„) - K0c + ^ \/l - r2sin^(0„ + (pc) - (2) 



- -yl- r2sin2(0„) . 

A solution of the function G{(j)c) for 0c ^ (0, 27r] corresponds to a collision of the particle 
with the moving wall and it is obtained numerically. 

Let us now consider the case where the particle leaves the collision zone, i.e. case 
(ii). The corresponding expressions are V* = Vn, AT„ = (pr + <Pc where 0^ corresponds 
to the elapsed time the particle spends travelling from the last hit with the moving wall, 
up to suffering an elastic reflection with the static wall and be reflected backwards, 
therefore until the entrance of the moving wall. Thus, 0^ is given by 



2 + U- Wl-r2sin2(0„) -ecos(0„)-e 
= ^ ^ . (3) 

The term 0c is numerically obtained from F{(j)c) = for 0c G [0, 2tx) where the function 
F[(j)c) is given by 



F(0c) = ecos(0„+0T+0c) + ^Vl - r^sm\(Pn + 0t + 0c)-^-e+K0c -(4) 

After some straightforward algebra, it is easy to show that the mapping ([1]) preserves 
the following phase space measure 



diJ, 



y + esin(0) I 1+ 



dVd^ . (5) 



'l — sin^(0) _ 

We stress that in the limit of r — 0, the results for the one-dimensional Fermi accelerator 
model are all recovered O |19] . 

Figure [2] shows the corresponding phase space obtained via iteration of mapping 
([T]) for the control parameter e = 0.01 and: (a) r = 0.1, (b) r = 0.3, (c) r = 0.6 and (d) 
r = 0.9. We can clearly see that the shape of the phase space changes as the control 
parameter r varies. On the other hand, the mixed form is preserved in the sense that a 
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Figure 2. Phase space generated from mapping ([T]) and control parameter e = 0.01 
and: (a) r = 0.1, (b) r = 0.3, (c) r = 0.6 and (d) r = 0.9. 



large chaotic sea, which surrounds KAM islands, is limited by a set of invariant spanning 
curves. One can also note that the position of the lowest invariant spanning curve raises 
as the control parameter r increases. In our simulations we will consider values for 
r that may approaches the unity, moreover in the range of r G [0,1]. However for a 
real experimental system, in which damping forces can not be neglected, such values 
for r have no much interest. This is mainly because the damping force can acquires 
larger values as compared to the component of the force with respect to the motion (for 
instance, it happens for large values of the angle ip) and therefore, lead the system to 
reach the rest. 
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Figure 3. Positive Lyapimov exponent for 10 different initial conditions randomly 
chosen along the chaotic sea in the low energy regime. The control parameters used 
were e = 0.01 and r = 0.1. 

The two natural questions that we are interested in concern on the properties of the 
chaotic sea (see Fig. [2]), hke the positive Lyapunov exponent and the average velocity 
of the particle. Our main goal is to describe their behaviour as function of the control 
parameter r. We think this study is of interest because the control parameter r directly 
controls the straight of a nonlinearity of the model. For small values of r, results of 
the FUM should be obtained. Moreover, we expect that the results obtained for r — ^ 1 
contribute towards a better understanding of this model for such a range of r and, 
in particular as we will see in Sec. [3l in a description of the phenomenon of Fermi 
acceleration. The behaviour of the Lyapunov exponent is described in this section while 
the average velocity is discussed in section [3l 

We concentrate to investigate the behaviour of the positive Lyapunov exponent for 
the chaotic sea. It is well known that the Lyapunov exponent is commonly used as 
a tool to characterise sensitivity to initial conditions. Figure [3] shows the asymptotic 
convergence of the positive Lyapunov exponent for the control parameters e = 0.01 
and r = 0.1. The ensemble average of 10 different initial conditions randomly chosen 
along the chaotic sea gives A = 0.82 ± 0.01, where the error 0.01 denotes the standard 
deviation of the ten samples. Each initial condition was iterated up to 10'' collisions 
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Figure 4. Positive Lyapimov exponent as function of r for a fixed e = 0.01. 



with the moving wall. The method used to obtain the Lyapunov exponents is described 
in the Appendix. 

Let us discuss the behaviour of the positive Lyapunov exponent as function of r. 
FigureHlshows the behaviour of A x r. We can see that, in the limit of r — > 0, the positive 
Lyapunov exponent recovers the value of the one-dimensional Fermi accelerator model 
[5|,[T9]. The positive Lyapunov exponent then grows slightly, having a maximum value 
around r = 0.15 and then decreases almost monotonically until around r = 0.5. Then 
it starts grow again until r ^ 0.7 when it suddenly decreases. Other abrupt change is 
observed for r ^ 0.89. Thus, the two main questions that arise from Fig. Hlare: (i) why 
does the positive Lyapunov exponent decreases in the average, instead of growth, as r 
raises? (ii) what is the explanation of the abrupt changes in A? 

The answer for question (i) comes from the shape of the function that describes 
the motion of the moving wall. Figure [5] shows four different plots of S{(j)) = 
rcos(0) + \Jl — r"^ sin^(0) (for sake of clarity, we show two periods in for S{(j))), 
which describes the motion of the moving wall for four different values of r namely: (a) 
r = 0.1, (b) r = 0.4, (c) r = 0.8 and (d) r = 0.999. For r = 0.1, the function looks like 
a cosine function but as r increases, the shape changes substantially. It thus become to 
have two regimes of variation where one of them is characterised by a constant plateau 
in the limit of r — 1, as can be seen in the range of G (7r/2, 37r/2) and a regime of fast 
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Figure 5. Plot of 5(0) for two perfods in </> considering the control parameters e = 0.01 
and: (a) r = 0.1, (b) r = 0.4, (c) r = 0.8 and (d) r = 0.999. 

variation, which is given by the complementary values of 0. In the limit case oir — 1, 
the function has indeed discontinuities in its derivative for two values of 0, namely 
= 7r/2 and = 37r/2. As we will see in the next section, the discontinuities for 
yield a profound consequence in the dynamics of the system thus leading the particle 
to exhibit unlimited energy growth. The large plateaus imply that the particle, when 
suffers collisions with the moving wall does not change substantially its velocity value 
since such plateaus lead in an "almost" null velocity for the moving wall. Thus the time 
that the particle spends until a next hit with the moving wall is almost the same as it 
spent in the previous collisions. It then implies that the particle, in a chaotic orbit, can 
experience many more collisions with the moving wall without substantially changing 
its energy as it would experiences if the plateaus were absent. Then, the form of 
for large values of r yields in reducing the chaoticity of the system in the average. 

We will now discuss a possible answer for question (ii), i.e., the explanation of the 
sudden changes in A. They are basically related to the destruction of the lowest energy 
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Figure 6. (a) Illustration of chaotic behaviour generated by 5 different initial 
conditions for the control parameters e = 0.01 and r = 0.69. (b) Iteration of a single 
initial condition, evidencing the chaotic behaviour, for e = 0.01 and r — 0.72. (c) 
Zoom in of part (a) where it is easy to see an invariant spanning curve and chaotic 
layers. 

invariant spanning curve together with a destruction of small chaotic layers. Considering 
the case of low energy and for r = 0.69, the system has a large chaotic sea (the black 
region of Fig. EJ^c)) which shares boundary with a thin chaotic layer (the red region 
of Fig. El^c)), as can be seen in Fig. [61 Above this chaotic layer, there is an invariant 
spanning curve (brown curve of Fig. [6](c)). Above yet of such curve we can see two thin 
chaotic layers (light blue and dark blue) and a relatively large chaotic sea (green region) 
surrounding a KAM island (see for instance Fig. El^c)). Each of these chaotic regions 
were characterised in terms of Lyapunov exponents. Their corresponding values were: 
for the black region we obtained Ab = 0.574(6); for the red region = 0.072(1); for the 
light blue Alb = 0.0110(4); for the dark blue Adb = 0.0152(2) and for the green region 
Ag = 0.218(1). For r = 0.72, which is quite close to r = 0.69, those regions shown in 
Fig. [6](a) were all merged into a single and large chaotic sea characterised by a positive 
Lyapunov exponent A = 0.46(5). The merged regions are shown in Fig. EJ^b). After 
the destruction of the thin structures shown in Fig. [6](a), the chaotic sea in the low 
energy region can spreads over a larger accessible region of the phase space. So we can 
consider that, after the transition, the positive Lyapunov exponent could be obtained 
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by an average of the previous values for the corresponding chaotic regions therefore 
taking into account the fraction of area occupied individually by each region. To check 
whether this supposition is correct, we have obtained the fraction of each chaotic region 
previous to the control parameter variation i.e. for r = 0.69. We have defined a grid of 
84 initial conditions for the 0-axis, limited to the interval G [0,27r), and 127 for the 
\^-axis considering the interval V e [0,0.38813]. The value 0.38813 corresponds to the 
higher value of the velocity obtained for the chaotic sea shown in the green region of Fig. 
MyC). Thus in the total, we considered 10668 different initial conditions. Each of them 
were evolved in time for n = 10'' collisions with the moving wall and their Lyapunov 
exponents evaluated. The obtained value was compared to the Lyapunov exponent of 
the chaotic regions so that we were able to compare the corresponding fraction of initial 
conditions which belongs to one region or to another one. Applying this procedure for all 
the chaotic regions shown in Fig. [6](c), we found that the black chaotic sea fills a fraction 
of Pb = 0.8261 of the entire chaotic region. The red region corresponds to = 0.0279, 
when the light blue has a fraction of Pib = 0.0103, the dark blue is Pdb = 0.0038 and 
finally the green region corresponds to a fraction of Pg = 0.1319. After the transition, 
we could assume that the positive Lyapunov exponent is obtained by 

A = AbPb + ArPr + AlbPlb + AdbPdb + AgPg . (6) 

Evaluating Eq. we found A = 0.50 which is a rather acceptable value as compared 
to the value obtained via numerical simulation of the chaotic sea after the destruction 
A = 0.46(5). 

The abrupt change in A around the value of r ^ 0.89 is explained by using the same 
arguments. We stress that similar results were observed in a rather distinct model |20j . 

3. A simplified version of the model and Fermi acceleration 

Other important conclusion that arises from the shape of the function 5(0) is related 
to the energy gain of the bouncing particle. As r approaches the unity, the variation of 
the moving wall position becomes more fast for specific ranges of 0. It implies that the 
particle can acquires, for specific ranges of 0, large values of velocity upon collisions with 
the moving wall for those regions of fast variation of 5'(0). Such a result can be seen 
in Fig. [2] by the position of the lowest energy invariant spanning curve which assume 
higher values as the control parameter r raises. To illustrate such an argument more 
clearly, it is important to look at the behaviour of the average velocity for sufficiently 
long time. To do this, we will make use of a simplification in the mapping ([T]) with 
the main goal of speeding up our numerical simulations and avoid solving the equations 
G(0c) = and P(0c) = 0. This simplification, which is commonly used in the literature 
(see Refs. [U El [TJ EH |22l ES]), consists in assume that both walls are fixed. However, 
when the particle hits one of them, it exchanges energy and momentum as if the wall 
were moving. This procedure retains the nonlinearity of the problem and yields a huge 
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advantage of avoid solving transcendental equations. The mapping is then given by 



T : 



0n+l 
K+1 



0„ + r^] mod(27r) 
Vn - 2esin(0„+i) 



^ _|_ rcos((^„+i) 



(7) 

Although this simplification brings the advantage of allowing very fast simulations as 
compared to those of the complete version, it also gives rise to a problem that we need 
to avoid. In the complete model, depending on the combination of both velocity and 
phase of the moving wall, it is possible for the particle, after suffering a collision with 
the time varying wall, to suffer a second successive collision before exiting the collision 
area, as well as possibly having a negative velocity following the first such a collision. In 
the simplified model, non-positive velocities are forbidden because they are equivalent 
to the particle travelling beyond the wall. In order to avoid such problems, if after the 
collision the particle has a negative velocity, we inject it back with the same modulus of 
velocity. Such a procedure is effected perfectly by the use of a modulus function. Note 
that the velocity of the particle is reversed by the modulus function only if, after the 
collision, the particle remains travelling in the negative direction. The modulus function 
has no effect on the motion of the particle if it moves in the positive direction after the 
collision. We stress that this approximation is valid only for small values of e. 

We now discuss the procedure used to obtain the average velocity of the particle. 
Firstly we obtain the average velocity proceeding with an average over the orbit, i.e. 
1 

V^=-I: V^, , (8) 

where j denotes the number of collisions. The second average is made in an ensemble 
of M = 10^ different initial conditions so that the average value is defined as 

F = V V^, , (9) 

with i corresponding to a sample in an ensemble of M different initial conditions. 

Figure [7](a) shows the behaviour of x n for the control parameters e = 0.01 and 
r = 0.01 considering n = 10^ iterations. We can see that the velocity of the particle 
grows for short iterations and then suddenly it bends towards a regime of saturation for 
sufficiently long time. Moreover, we are interested in the behaviour of the asymptotic 
values of V. Figure [7](b) shows the behaviour of V for long iterations, which we will 
referr to it as V^at- The behaviour of Kat can be described for two different ranges of r. 
The first range is for r < 0.9 where we can see that V^at is almost constant for the region 
of r G [0,0.9). The second range of r is considered for r > 0.9. It is important to say 
that, when r = 1 the expression of the "moving wall" velocity presents discontinuities 
for (j) = n/2 and (p = 3tt/2. Thus, it is convenient to define a new parameter fi = 1 — r 
and then study the behaviour of V^at as a function of fi. This new control parameter has 
a practical application because it brings the criticality of the model to /i = 0. Figure [8] 
shows the behaviour of V^at x for a fixed control parameter e = 0.01. We can describe 
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Figure 7. (a) Behaviour oi V x n for tlie control parameters r = 0.01 and e = 0.01. 
(b) Plot of Fsat X r for a fixed e = 0.01. 

such a behaviour according to 

^sat CX /i" . (10) 

After fitting a power law in Fig. [HI we obtain a = —0.460(6) = —0.5. This kind 
of behaviour for Vsa,t confirms that, in the limit of r 1, the present model shows 
the phenomenon of Fermi acceleration. This result has a clear explanation in terms 
of the KAM theorem. As it was discussed by Lichtenberg, Lieberman and Cohen [4], 
if the expression of the periodic wall's velocity has a sufficient number of continuous 
derivatives, then it is possible to obtain invariant spanning curves separating different 
portions of the phase space. Particularly they are useful to prevent unlimited energy 
growth for a bouncing particle. Moreover, it was estimated by Moser |23] that three 
continuous derivatives for the moving wall velocity is a sufficient condition for the 
existence of KAM surfaces (by instance, invariant spanning curves). However, in the 
present model and considering the limit of r — > 1, the expression of the velocity shows 
discontinuities for (p = tt /2 and = 37r/2, so that the invariant spanning curves are not 
observed for r = 1 and consequently, Fermi acceleration is present. 

4. Conclusions 
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In summary, we have studied a non dissipative version of a classical bouncing ball 
model under the presence of two nonlinearities. Our results show that, as one of the two 
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Figure 8. Plot of V^at x n for a fixed e = 0.01. A power law fitting gives that 
a = -0.460(6) ^ -0.5. 

control parameters varies, the positive Lyapunov exponent diminish in the average and 
experiences sudden changes. We have explained such a behaviour by the shape of the 
moving wall and due to the destruction of invariant spanning curves and thin chaotic 
layers. We have also shown that in the limit of r 1, the present model exhibits 
unlimited energy growth. This phenomenon was explained by using a discontinuity of 
the moving wall's velocity. 
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Appendix 

In this section, we briefly discuss the procedure used to obtain the Lyapunov exponents. 
In effect, the procedure consists in evolving the system over a long time from two slightly 
different initial conditions. If the two trajectories diverge exponentially in time the orbit 
is called chaotic, and the Lyapunov exponent obtained is positive. If the Lyapunov 
exponent is negative, the orbit may be either periodic or quasi-periodic. Let us now 
describe the procedure used to obtain the Lyapunov exponents numerically. They are 
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defined as ^ ^ 

A,- = lim In I A,- 1, ?' = 1,2, 

where Aj are the eigenvalues of M = 11^=1 <^fc(^A:5 ^fc) and Jk is the Jacobian matrix 
evaluated over the orbit (Vfc,0fc). The Jacobian matrix is defined as 



J 



dVn + l dVn^ 

dVn d4>n 

d4>n+l d4>„-{ 

dV„ (9</<„ 



In order to evaluate the eigenvalues of M, we use the fact that J can be written as a 
product of J = BT, where is an orthogonal matrix and T is a right upper triangular 
one. We now define the elements of these matrices as 



e 



cos(6') — sin(6') 
sin(^^) cos(6') 







Tl2 \ 




V 


T22 1 



Since M is defined as M = JnJn-i ■ ■ ■ J2J1, we can introduce the identity operator, 
rewrite M as M = J„ J„_i ... 72616]'^ Ji, and define Qi^Ji = Ti. The product 
J2Q1 defines a new matrix J|. In a following step, we may write M as M = 
JnJn-i ■ ■ ■ <^30202 ^<^2T'i- The Same procedure yields T2 = 02^^J2- The problem is 
thus reduced to the evaluation of the diagonal elements of Tj : T^^, T22. Using the 9 and 
T matrices, we find the eigenvalues of M, given by 

^ _ ill + J2I _ JIIJ22 - jl2j21 

~ ; ^ =; -'22 — / ^ ^ — - 

ill + i21 Y ill + i21 

We can then evaluate the Lyapunov exponent using the relation 

n 1 

A, = lim 5:-ln|T,^.|, J = 1,2. 

k=l 

It is interesting to observe that Ai = — A2, because the map is measure-preserving. 
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